"""
    AdvancedSolvers

Advanced optimization algorithms module for GSI Core Analysis.
This module provides sophisticated iterative solvers ported from GSI Fortran:

# Solver Categories

## 1. Lanczos Methods (`LanczosSolver`)
- Complete Lanczos tridiagonalization process
- Implicit restart mechanisms for large-scale problems
- Ritz value computation and selection
- Memory-efficient vector storage with reorthogonalization
- Preconditioning support for enhanced convergence

## 2. BiCG-Lanczos Methods (`BiCGLanczos`) 
- BiCG-Lanczos for non-symmetric systems
- Bi-orthogonal Lanczos process with left/right vectors
- Multiple algorithmic variants:
  - BiCGStab: Bi-Conjugate Gradient Stabilized
  - BiCGStab(l): Higher-order stabilization
  - QMR: Quasi-Minimal Residual approach
  - Look-ahead BiCG: Breakdown-free variant
- Robust error handling and breakdown recovery

## 3. Quasi-Newton Methods (`QuasiNewton`)
- L-BFGS: Limited-memory BFGS with variable memory
- Multiple line search strategies (Armijo, Wolfe, More-Thuente)
- Trust region methods for global convergence
- Hessian approximation with efficient storage
- Preconditioning and scaling strategies

# Mathematical Framework

These solvers address the core GSI variational analysis equation:
```
minimize J(x) = ½(x-xᵦ)ᵀB⁻¹(x-xᵦ) + ½(H(x)-y)ᵀR⁻¹(H(x)-y)
```

Where the solvers handle the iterative minimization through:
- **Lanczos**: Eigenvalue problems and preconditioning
- **BiCG-Lanczos**: Non-symmetric linear systems  
- **Quasi-Newton**: Nonlinear optimization with approximate Hessian

# Usage Examples

## Lanczos Solver
```julia
using GSICoreAnalysis.AdvancedSolvers

config = LanczosConfig(max_iter=50, tolerance=1e-8, precondition=true)
result = lanczos_solve!(x, A, b, config)
```

## BiCG-Lanczos Solver  
```julia
config = BiCGLanczosConfig(max_iter=100, variant=BiCG_STAB, tolerance=1e-10)
result = bicg_lanczos_solve!(x, A, b, config)
```

## Quasi-Newton Solver
```julia
config = QuasiNewtonConfig(max_iter=50, memory_size=10, line_search_type=WOLFE)
result = lbfgs_solve!(x, objective, gradient!, config)
```

# Integration with GSI Framework

These solvers integrate seamlessly with:
- `ControlVectors`: For state vector operations
- `CostFunctions`: For objective function evaluation
- `BackgroundError`: For preconditioning operations
- `StateVectors`: For efficient memory management

# Performance Considerations

- **Memory efficiency**: Limited storage for large-scale problems
- **Numerical stability**: Robust algorithms with breakdown recovery
- **Convergence**: Advanced stopping criteria and adaptive methods
- **Parallelization**: Support for distributed computing environments

The implementations maintain mathematical equivalence with GSI Fortran
while leveraging Julia's performance and expressiveness.
"""
module AdvancedSolvers

using LinearAlgebra
using SparseArrays
using Printf

# Include submodules
include("LanczosSolver.jl")
include("BiCGLanczos.jl") 
include("QuasiNewton.jl")

# Import and re-export from submodules
using .LanczosSolver
using .BiCGLanczos
using .QuasiNewton

# Export Lanczos solver
export LanczosConfig, LanczosResult, LanczosState
export lanczos_solve!, setup_lanczos_precond, apply_lanczos_precond
export compute_ritz_values, lanczos_convergence_check

# Export BiCG-Lanczos solver
export BiCGLanczosConfig, BiCGLanczosResult, BiCGLanczosState
export bicg_lanczos_solve!, bicgstab_solve!, qmr_solve!
export setup_bicg_precond, apply_bicg_precond
export BiCGVariant

# Export Quasi-Newton solver
export QuasiNewtonConfig, QuasiNewtonResult, QuasiNewtonState
export lbfgs_solve!, bfgs_solve!, dfp_solve!
export LineSearchType, TrustRegionType
export setup_quasi_newton, apply_quasi_newton_precond

# Export enumerations
export BiCG_STANDARD, BiCG_STAB, BiCG_STAB_L, QMR, LOOK_AHEAD_BICG
export ARMIJO, WOLFE, STRONG_WOLFE, MORE_THUENTE, BACKTRACKING
export NO_TRUST_REGION, CAUCHY_POINT, DOGLEG, STEIHAUG_CG

# Convenience functions for unified interface
"""
    create_default_configs(problem_size::Int, precision::Type{T}=Float64) where T

Create default configurations for all three solver types.

Returns a named tuple with optimized configurations for typical GSI problems.
"""
function create_default_configs(problem_size::Int, precision::Type{T}=Float64) where T
    # Estimate appropriate parameters based on problem size
    memory_size = min(20, max(5, problem_size ÷ 100))
    max_iter = min(200, max(50, problem_size ÷ 10))
    
    lanczos_config = LanczosConfig{T}(
        max_iter = max_iter,
        tolerance = T(1e-8),
        precondition = true,
        max_precond_vecs = memory_size,
        verbose = false
    )
    
    bicg_config = BiCGLanczosConfig{T}(
        max_iter = max_iter,
        tolerance = T(1e-8),
        variant = BiCG_STAB,
        precondition = true,
        verbose = false
    )
    
    qnewton_config = QuasiNewtonConfig{T}(
        max_iter = max_iter ÷ 2,  # Fewer iterations for nonlinear
        tolerance_grad = T(1e-6),
        memory_size = memory_size,
        line_search_type = WOLFE,
        verbose = false
    )
    
    return (
        lanczos = lanczos_config,
        bicg = bicg_config,
        quasi_newton = qnewton_config
    )
end

"""
    solve_advanced!(x, problem, method::Symbol, config=nothing; kwargs...)

Unified interface for all advanced solvers.

# Arguments
- `x`: Solution vector (modified in place)
- `problem`: Problem definition (matrix, operator, or objective function)
- `method`: Solver method (:lanczos, :bicg_stab, :lbfgs)
- `config`: Solver configuration (optional, defaults created if not provided)

# Examples
```julia
# Linear system with Lanczos
result = solve_advanced!(x, A, :lanczos, config; b=b)

# Linear system with BiCGStab
result = solve_advanced!(x, A, :bicg_stab, config; b=b)

# Nonlinear optimization with L-BFGS
result = solve_advanced!(x, obj_func, :lbfgs, config; gradient!=grad_func!)
```
"""
function solve_advanced!(x::VT, problem, method::Symbol, config=nothing; kwargs...) where VT
    
    n = length(x)
    T = eltype(x)
    
    if method == :lanczos
        # Lanczos method for symmetric eigenvalue problems or linear systems
        if config === nothing
            config = LanczosConfig{T}(max_iter=min(100, n), tolerance=T(1e-8))
        end
        
        b = get(kwargs, :b, nothing)
        if b === nothing
            error("Lanczos solver requires right-hand side vector 'b'")
        end
        
        return lanczos_solve!(x, problem, b, config; kwargs...)
        
    elseif method == :bicg_stab
        # BiCGStab for non-symmetric linear systems
        if config === nothing
            config = BiCGLanczosConfig{T}(max_iter=min(100, n), variant=BiCG_STAB, tolerance=T(1e-8))
        end
        
        b = get(kwargs, :b, nothing)
        if b === nothing
            error("BiCG solver requires right-hand side vector 'b'")
        end
        
        return bicgstab_solve!(x, problem, b, config; kwargs...)
        
    elseif method == :lbfgs
        # L-BFGS for nonlinear optimization
        if config === nothing
            config = QuasiNewtonConfig{T}(max_iter=min(50, n÷2), tolerance_grad=T(1e-6))
        end
        
        gradient! = get(kwargs, :gradient!, nothing)
        if gradient! === nothing
            error("L-BFGS solver requires gradient function 'gradient!'")
        end
        
        return lbfgs_solve!(x, problem, gradient!, config; kwargs...)
        
    else
        error("Unknown solver method: $method. Available: :lanczos, :bicg_stab, :lbfgs")
    end
end

"""
    benchmark_solvers(problem_size::Int, num_trials::Int=5)

Benchmark all three solver types on synthetic problems.
"""
function benchmark_solvers(problem_size::Int, num_trials::Int=5)
    println("=== Advanced Solvers Benchmark ===")
    println("Problem size: $problem_size")
    println("Number of trials: $num_trials")
    println()
    
    # Generate test problem
    A = randn(problem_size, problem_size)
    A = A + A'  # Make symmetric for Lanczos
    b = randn(problem_size)
    x_true = A \ b
    
    results = Dict()
    
    # Benchmark Lanczos
    println("Testing Lanczos solver...")
    times_lanczos = Float64[]
    errors_lanczos = Float64[]
    
    config_lanczos = LanczosConfig(max_iter=min(50, problem_size), tolerance=1e-10)
    
    for trial = 1:num_trials
        x = zeros(problem_size)
        start_time = time()
        result = lanczos_solve!(x, A, b, config_lanczos)
        elapsed = time() - start_time
        
        push!(times_lanczos, elapsed)
        push!(errors_lanczos, norm(x - x_true))
    end
    
    results[:lanczos] = (
        mean_time = sum(times_lanczos) / num_trials,
        mean_error = sum(errors_lanczos) / num_trials,
        times = times_lanczos,
        errors = errors_lanczos
    )
    
    # Benchmark BiCGStab
    println("Testing BiCGStab solver...")
    times_bicg = Float64[]
    errors_bicg = Float64[]
    
    config_bicg = BiCGLanczosConfig(max_iter=min(100, problem_size), variant=BiCG_STAB, tolerance=1e-10)
    
    for trial = 1:num_trials
        x = zeros(problem_size)
        start_time = time()
        result = bicgstab_solve!(x, A, b, config_bicg)
        elapsed = time() - start_time
        
        push!(times_bicg, elapsed)
        push!(errors_bicg, norm(x - x_true))
    end
    
    results[:bicg_stab] = (
        mean_time = sum(times_bicg) / num_trials,
        mean_error = sum(errors_bicg) / num_trials,
        times = times_bicg,
        errors = errors_bicg
    )
    
    # Print summary
    println("\n=== Benchmark Results ===")
    for (method, data) in results
        @printf("%-12s: %.4f ± %.4f sec, error: %.2e ± %.2e\n", 
                string(method), 
                data.mean_time, 
                std(data.times),
                data.mean_error, 
                std(data.errors))
    end
    
    return results
end

# Print module summary when loaded
function __init__()
    println("AdvancedSolvers module loaded successfully")
    println("Available solvers: Lanczos, BiCG-Lanczos, Quasi-Newton")
    println("Use create_default_configs() for optimized configurations")
end

end # module AdvancedSolvers